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Abstract — Information about the strength of gas sources in 
buildings has a number of applications in the area of building 
automation and control, including temperature and ventilation 
control, fire detection and security systems. Here, we consider the 
problem of estimating the strength of a gas source in an enclosure 
when some of the parameters of the gas transport process are 
unknown. Traditionally, these problems are either solved by the 
Maximum-Likelihood (ML) method which is accurate but com- 
putationally intense, or by Recursive Least Squares (RLS, also 
Kalman) filtering which is simpler but less accurate. In this paper, 
we suggest a different statistical estimation procedure based on 
the concept of Method of Moments. We outline techniques that 
make this procedure computationally efficient and amenable for 
recursive implementation. We provide a comparative analysis of 
our proposed method based on experimental results as well as 
Monte-Carlo simulations. When used with the building control 
systems, these algorithms can estimate the gaseous strength in a 
room both quickly and accurately, and can potentially provide 
improved indoor air quality in an efficient manner. 

Index Terms — Occupancy Estimation, Parameter Estimation, 
Non-linear Regression, Monomolecular Growth Curve, Method 
of Moments 



I. Introduction 

Recent studies have indicated that indoor air quality (IAQ) 
has a significant effect on the productivity and health of the 
occupants, e.g., for the U.S., the estimated potential annual 
savings and productivity gains due to improved indoor air 
quality are $30 billion to $170 billion [1]. Consequently 
a plethora of research is being aimed at improving indoor 
environmental systems, yet making sure that they are en- 
ergy efficient and sustainable [2|. One such energy-efficient 
technique is demand controlled ventilation (DCV). DCV is 
a well documented method that reduces heating and cooling 
requirements in buildings by adjusting ventilation rates in 
response to occupancy 0, 0. Determination of occupancy 
inside a space is essential to implement DCV, and this can 
be achieved by processing data from various sources like gas 
sensors, cameras, microphones, or may be even via manual 
logbooks and electronic calenders. Carbon dioxide (CO2 ) 
sensors are widely used because of cost considerations and 
also due to privacy concerns. The basis of using CO2 for occu- 
pancy estimation is established in well-quantified principles of 
human physiology. All humans, given a similar activity level, 
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exhale CO2 at a predictable rate based on occupant age and 
activity level, e.g., 0, |d|. 

The mass-balance model for CO2 generation in buildings 
is well established, e.g., 0, 0. An indoor CO2 measure- 
ment provides a dynamic measure of the balance between 
CO2 generation in the space, representing occupancy and the 
amount of low CO2 concentration in outside air introduced 
for ventilation (also known as inflow-rate). As a result, the 
indoor CO2 increases gradually until it achieves an equilibrium 
point where the CO2 produced by people is in balance with 
the dilution rate. A DCV based control strategy can either 
act on the transient (pre-equilibrium) data or wait for the 
steady state to be reached |9| before estimating the occupancy. 
Since at low inflow-rates it takes longer to reach equilibrium, 
steady-state algorithms usually involve significant delay 0. In 
this paper, we consider using the transient data for occupancy 
estimation. 

The problem of estimating the CO2 generation rate (or 
occupancy) and/or inflow parameters using transient data has 
been considered by several researchers. Most of the available 
techniques use the differential model of the mass-balance 
equations for estimation followed by a filtering operation to 
smooth the estimates. For example, an autoregressive filter was 
used by [3| and Kalman filters were used by [8| and |10|. 
While these techniques are very effective when dealing with 
dynamic situations, they often produce sub-optimal estimates 
when the occupancy is constant or slowly varying, as would 
be the case for meetings, classrooms, etc. As noted in 0, 
this deterioration of performance is due to the use of time- 
derivatives of concentrations. Since the derivatives cannot 
be measured directly, they must be approximated, and this 
approximation amplifies high frequency noise. 

When the occupancy in a space is constant, the integral 
model of the mass-balance equations resembles a simple 
growth curve, whose parameters are related to the occupancy 
and inflow-rate. Such a model is known to be adequately 
accurate for several applications 0. In this paper, we estimate 
the occupancy and inflow-rate by estimating the parameters 
of the growth curve and propose a new approach based on 
Method of Moments (MME) HH. Since our MME approach 
involves the use of time integrals which are more robust 
to noise, it performs better in terms of estimation accuracy. 
This fact is demonstrated using both theoretical analysis and 
experimental data. We also examine the performance of our 
algorithm in cases when the occupancy is slowly varying or 
when there is uncertainty regarding the activity levels of the 
occupants. A preliminary version of this work was presented 

Though many researchers have suggested procedures to esti- 
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mate the CO2 generation rate, performance analysis regarding 
the uncertainty of estimates has been scarce. One principal 
goal of this paper is to derive performance bounds that will be 
useful to practitioners. It is well known in statistical estimation 
literature (e.g., ifTTI . Ifl3ll ) that maximum likelihood estimator 
(MLE) has the minimum expected variance and is asymp- 
totically optimal, i.e., attains the Cramer-Rao Lower bound 
(CRLB). In this paper, we derive approximate closed form 
expressions for the CRLB and also corresponding performance 
metrics for the KF based procedure (presented in J8J) and the 
Method of Moments based approach. We demonstrate that for 
the transient region of operation, the MME estimator is more 
accurate than the KF counterpart. This improved accuracy can 
translate to significant reductions in response times and/or 
reduction in number of sensors. 

The rest of the paper is organized as follows. In Section 
|ll] we describe the model of the dynamics inside a well- 
mixed space, formulate our estimation problem, and introduce 
certain approximations that will aid our theoretical analysis. In 



Section III we discuss the MLE and Kalman Filtering based 



techniques and derive the CRLB and KF performance metrics. 



In Section IV we describe our time-integral based approach, 
comment on computational issues and derive the expected 
variance. In Section [V] we describe an experimental setup 
that assumes near ideal conditions and compare the estimation 
performances. We also present performance comparisons using 
theoretical results and Monte-Carlo simulations. Concluding 
remarks are provided in Section fVl] 



II. Problem Formulation 

The model of the dynamics inside a well-mixed space 
is described in this section. This is identical to the well- 
known single-zone model (e.g., |7|,|8|,|14|) and we make the 
following assumptions - the mass of the air in the space is 
constant, and the concentration distribution in the space is 
spatially uniform. Consider the schematic diagram of an indoor 
space in Figure [T] 



Q, C c 



v 



C t 



M 




Fig. 1. Schematic diagram of an ideal indoor space. N is the number of 
persons that emit CO2 at a total rate of cN kg/sec. Ct,Co are the CO2 
concentration (ppm) of indoor and supply air. 



Let N denote the number of occupants, c (kg/sec) the CO2 
generation rate per person, M (kg) the mass of air inside the 
room and t (sec) a particular time instant. Let Q (kg/sec), 
Ct (ppm or mg/kg) and Co (ppm) be the airflow rate and 
CO2 concentration of indoor and supply air respectively. The 
number of occupants N and the airflow rate Q are assumed 
to be unknown quantities that need to be estimated. 

A macroscopic mass balance in such a space results in the 
following mass continuity equation that is widely used by 
both the academic (8), iTPfl and industrial |[T5l communities 
to describe gas transport processes in an indoor environment 



MC t = -Q(C t - C ) + cN. 



(1) 



Assuming that there were no gas sources in the room already, 
the integration yields 



C t = C + (cN/Q) (1 - exp(-iQ/M)) . 



(2) 



Now we reduce the number of unknowns in the model. We 
note here that c and M are constants and hence can be either 
known beforehand or obtained via training. We assume Co to 
be constant over the interval of interest. Measuring the supply 
duct concentration, we can obtain a fairly accurate estimate 
of Co and we treat it as known. Let the sampling interval be 
T s and the total duration of measurements be T. This means 
that the i th sample is obtained at time t = iT s , where i = 
1, 2, ... n and T = nT s . Defining a t = C t — Co, our models 
in differential and integral domains are 

den 
dt 

ai {N, Q) = (cN/Q) (1 - exp(-iT s Q/M)) , (4) 



cN — Qai , and 



(3) 



respectively. We will use the notations <Zj and a,-(iV, Q) 
interchangeably throughout this paper. In expressions where 
the dependence on N and Q needs to be stressed explicitly, we 
use the longer version of the notation. Next, we characterize 
the measurement errors. While obtaining the measurement for 
C we assume that the CO2 sensor is corrupted by a time- 
independent additive zero-mean Gaussian noise with vari- 
ance a 2 . We assume a 2 to be known since the measurement 
error of sensor noise can be learnt beforehand. Let yi denote 
the noise-added observation for a,, i.e., yi — at +£,-. In vector 
notations, let a = [ai, a 2 , . . ■ , a„]', y = [y 1 ,y 2 ,--- y n ]' and 
e = [ei, £2, • • ■ , e„]', so that 



y 



(5) 



Our goal here is to estimate 6 := [N, Q]', which includes the 
estimated number of occupants N. 

The time-series of the indoor CO2 concentration, as pre- 
dicted by the mass balance equation stated above, resembles 
a growth curve that increases gradually before settling down 
to an equilibrium. For example, the indoor CO2 concentration 
during some of the meeting times are showrQ in Figure [2] 
In all cases, there is a region before equilibrium where the 
measurements increase steadily. We refer to this part as the 
transient region. 

'This dataset was acquired during a collaborative research project involving 
United Technologies Research Center and Syracuse University. Concentration 
values were sampled every minute from commercially available TI-4GS- 
22C06 sensors with rated accuracy of ±75 ppm. 
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Fig. 3. Growth curve predicted by the ideal indoor space model (solid line) 
Fig. 2. Time-response of the indoor CO2 concentration during some of the and polynomial approximations of increasingly higher order (dashed lines), 
meeting times. 



For the purpose of timely ventilation, the controller needs 
an estimate of N as early as possible, preferably before 
equilibrium is reached. In other words, our desired operating 
region is the transient region of the curve. We note that the 
transient region can be approximated as a p th order polynomial, 



(cN/Q) 1 - 



E 

,3=0 



(-iTsQ/Aiy 



(6) 



We define a parameter that is indicative of the transient region, 
K, as follows, 



K 



QT 

M ' 



(7) 



which is the total amount of fresh air injected normalized 
with respect to the volume of the room. Figure [3] gives us 
an idea on the value of p that describes with reasonable 
accuracy. We have considered the parameters corresponding to 
the experimental setup described later in V-A i.e., one person 
in a 780 cu.ft. room ventilated by fresh air at 28 cu.ft./min. The 
resulting rise time is M/(QT S ) = 30 samples and equilibrium 
value cN/Q — 530 ppm. 



III. Related Work 

As discussed in Section [II] our goal is to estimate the 
occupancy N and inflow-rate Q from the observation sequence 
y = a + e. In the non-linear regression literature (e.g., 
lfTD .fl6:h. the time series a given by Q is also known 
as Monomolecular Growth Curve. We discuss two existing 
approaches to solve this estimation problem and analytically 
evaluate their performance. 



A. Maximum Likelihood Estimation and CRLB 

In the MLE approach ifTTl . we minimize the log-likelihood 
function C for the measurements y, i.e., 

C(N,Q) = lnf(y;N,Q) (8) 

1 n 

= -nln(V2W) - ^^(^- ai (Ar,Q)) 2 , 

i—l 

(9) 

and obtain the estimates N and Q. Since are i.i.d. Gaussian, 
MLE is equivalent to minimizing the sum of squares 



(10) 



i=l 



The asymptotic variance of MLE attains the Cramer-Rao 
Lower Bound QT), which is the inverse of the Fisher infor- 
mation matrix 



T(N,Q) 



i /T(<h±\ 2 V daj dm \ 
J_ll^\dNJ Z^dNdQ] /in 

ON dQ 2^i\dN) / 



With the assumptions that (1) T is chosen such that K = 
QT /M is small (e.g., for our experimental setup, QT/M « 
0.5 corresponding to n — 50 samples) and (2) a second order 
p = 2 is sufficient to describe the signal in (|6j, we obtain the 
following approximations for the partial derivatives 



daJdN « (cT s /M)i, and 
da t /dQ^ {cN/2){T s /Mfi 2 . 



(12a) 
(12b) 



With another assumption that (3) the number of measurements 
n for yi is large, the CRLB can be derived as 



CRLB 



A8a 2 T s M 2 



(13) 
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When higher order polynomials are required to approximate 
the signal, using a similar approach, the CRLB can be more 
accurately expressed in the following form, 



CRLB 



48cr 2 T,M 2 



1 + aiK ■ 



K 1 



(14) 



where (1) K = QT/A1, (2) higher j implies more accuracy 
and (3) a>j are appropriate constants, e.g., the first three 
terms can be derived as a.\ — 2/3, a 2 = 377/945 and 
«3 = 2411/11340. We skip the derivations for the sake of 
brevity. The CRLB will be compared later with the perfor- 
mance metrics for other techniques in Section [V] 

B. Recursive Least Squares based Estimation 

The Kalman Filter based approach detailed in [8 1 simplifies 
to a Recursive Least Squares (RLS) algorithm (e.g., fPTl ') with 
the assumption of constant parameters. Based on the differen- 
tial model for signal in ([3]), the RLS model is predicated on 
the difference equation 

Vi ~ Vi-i = {T s /M)(cN - go*) + - e ( _i), (15) 

where e, represent the additive noise (observation error) and 
Ui represents the measurements of the room concentration. 
Equation ( fl"5j ) above represents a set of n equations for the 
Ti+l measurements for i = 0, 1, . . . , n. Denoting 7jj as 



Vi 



(16) 



the equations can be represented as Y = X0 + r), where 
= [N, Q]' and X is the matrix containing the signal terms 



X 



M 



c c 



c 

«„ 



(17) 



From the theory of Generalized Least Squares (e.g., lfTT1 ). we 
know that the LS estimate and the expected variance are given 
by 

= [X'V^Xy 1 X'V^Y, and (18) 
Var(0) = [X'V^xy 1 , (19) 

where V = Var(r7). We are interested in the variance of N, 
which is given by 



Var(iV) = [[X'y- 1 ^] 1 



(20) 



The matrices V and V~ x need to be described further. Note 
that the time correlation of r]i results in this n x n tridiagonal 
matrix 



V = a I 



-1 
2 



(21) 



The inverse of this tridiagonal matrix can be written in closed 
form. Bearing in mind that V~ x is symmetric, here we write 



only the upper diagonal terms 



a 2 {n + 1) 
n n — 1 
2(n- 1) 



I 

2i 



(n ~ i + l)i 



1 

2 

7i — i + 1 



(22) 



Let o\ LS denote the variance given by ( pO] ). Expressing the 
variance in a form similar to the CRLB ( fT4] > in Section |III-A| 
we obtain 



2 

a RLS 



192ct 2 T,A/ 2 



1 + PiK 



K 3 



where (1) K = QT/M, (2) higher j implies more accuracy 
and (3) j3j are appropriate constants, e.g., first three terms can 
be derived as ft = 3/8, /3 2 = 169/1120 and /3 3 = 57/1120. 
We skip the derivations for the sake of brevity. This variance 
can be compared with the CRLB ( fL4] > at this point. 

IV. Method of Moments based Estimation 

In this paper, we adopt a Method of Moments (e.g., 
lfTTIl . lfT7l ) based approach to estimate the occupancy and 
inflow-rate parameters. Unlike least squares approximation, 
the idea here is to match some theoretical and empirical 
properties of the model and solve the resulting equations to 
obtain the estimates. 

For m > 1, let us define by A m ^{N, Q) the theoretical 
integral of m power of at(N, Q) in the time duration [0, T], 



<N,Q) 



(a t (N, Q)) m dt. 



(24) 



It is possible to obtain an empirical integral by approximating 
Jo y m (t)dt with the help of the values {yo, j/j, . . . y n } at 
the set of time instants {0, T s , 2T S , . . . nT s }. We define the 
approximation as Y m ,T, 



Y m , T = I y m (t)dt, 



(25) 



so that the theoretical and empirical moments can be matched 



iN,Q) =Y„ 



(26) 



Since we have to estimate two parameters, we can obtain 
two equations by matching any two moments. We must note 
here that the choice of moment functions is not unique and one 
has to base the choice on three considerations - (1) existence 
of a unique solution, (2) computational complexity in solving 
the equations and (3) expected variance of estimation error. In 
Section |IV-A| we choose a family of moment functions that 
have a unique solution and the solution of which is amenable 
for recursive implementation as sample size increases. In 



Section IV-B we find the expected variance of the estimator. 
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A. Solution of Moment Equations 

We consider the system of two equations matching the 1 st 
and any m th (m > 2) moments. We use this to solve for our 
MME estimates N and Q. 

Ax, T {N, Q) = Y ltT , and A ljT (N, Q) = F m , T) m > 2 

(27) 

From our model description in Q, A„ u t(N, Q) can be 
described as 



f) T L <'--->"*■ (2., 

where r = QT/M, z = t/T. 



A m , T (N,Q) 



From the expressions for A m ^T in ( |28] i, one can observe that 
these are non-linear equations in N and Q. We demonstrate 
that they have a unique solution by presenting an approach in 
which these solutions can be obtained. We define a function 
G m (r), obtained by eliminating N from the expressions for 
A\t and A m t, as follows 



G m {r) 



(Ai, T ) m 



(29) 



It is easy to show that G m (r) is a decreasing convex function 
in (0,oo) for m — 2 and this property appears to be true for 
all G m (r) for m > 2. Monotonicity of G m (r) ensures the 
existence of an inverse while convexity of G m (r) ensures that 
the inverse can be computed quickly using any gradient-based 
numerical method. The solution for inflow-rate is 



Q 



M 



Tm-ly 
1 rr 



(30) 



We use the facts (1) F 1;T = Ai >T (N, Q) (the first moment 
condition) and (2) A 1>T (N,Q) = NA 1>T (l,Q) (follows from 
(|28|i), to obtain the estimate of the number of occupants 

M Fit 



N = 



(31) 



The existence of a unique solution means that our MME 
estimator is well defined. From the theory of Method of 
Moments [17|, this uniqueness property also implies statistical 
consistency of the estimators, i.e., given a large number of 
samples, the estimates converge to the true values. 

B. Expected variance of occupancy 

To compare the performance of the MME method with that 
of MLE or RLS, we derive the expected variance of our MME 
estimator for the number of occupants N . First, let us define 
the following vectors 

A = (A 1>T) A mjT )', and Y = (Y hT , Y m , T )' (32) 

so that our MME estimates satisfy 

A$) = Y, (33) 



where 6 = [N,Q]', For small estimation errors (i.e., small 
values of a 2 ), one can obtain the Taylor series approximation 
of A(0) around 6 



This can be rewritten in the following form 

-l 



6 6 

We denote 
H 



dA{6) 
86' 



A(6)-A(6] 



dA{6) 



86' 



and S = Var 



(A(£ 



(34) 



(35) 



(36) 



Asymptotic normality of A{6), which we show later, and ( [35] 
then implies Var(6>) ps HUH' . We are interested in Var(iV) 



which is actually 



Var(6>) 

Var(TV) 



and hence 



\HSH' 1 



(37) 



Next, we derive simplified expressions for S and H. 

We note from ( |33| l and ( f3"6] > that S = VarF. If the sampling 
interval is small enough compared to the rate of change of y t , 
we can use constant interpolants Y m> T = T s uT ■ Therefore, 
Var (F m , T ) = T s 2 ]T Var (j/f 1 ). Hence, for m = 1, 



Var(F 1 ,„)=T s 2 ^Vai-( 2/j )=T ; 



2 2 

no . 



(38) 



For calculating the variance for m > 2, we consider the 
expansion of y™ = (a, + £j) m and write their expectation 
as follows 



Vi 



a,- 



(39) 



m-2 2 

a„ cr 



The influence of the noise terms ej , € i , 



1 on the overall 
variance decays exponentially since they are multiplied by 
decaying factors a™ -1 , a™ -2 , etc. Hence, we can approximate 
y™ 1 by retaining only the first two terms from |39|, which 
means that y™ 1 is approximately Gaussian with mean a™ and 
variance 

Var(yn=Var( 2/ r-E(yr)) 



= Var (ma[ 
« Var fmol 



+ o(C 2 ) ( e ?-^ 2 ) + ---) 



Next, we can complete our derivation 



Var(F m , T ) =T 2 > Var (ma?' > , 



2m-2 



Using exactly similar arguments, the covariance term is 

Cov (Y hT , Y m>T ) w T 2 ma 2 V a™" 1 , 



(40) 



(41) 



(42) 



so that A{6) is approximately normal with covariance matrix 



(cNTs/M)™- 1 ^ 



(cNT s /M) r 



m 2 (cNT e /M) 2m ~ 2n 



2m-l 

(43) 
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with the assumptions that (1) the signal in |6]) is adequately 
described by p = 2 and (2) number of observations n is large. 

Assuming small QT/M, we obtain the following first order 
approximation of A m ^ 



A, 



M 



,771+1 



mQT s n m+2 
2M m + 2 



(44) 



from which we can obtain H. Using the expressions for H 
and X in ( [37] i, we obtain 



Var(TV) 



M 



4(m+ l) 2 (12 + m) cr 2 
m(2m — 1) n 3 



(45) 



For an MME algorithm with moment rn = 2, we denote the 



variance as a 



MME,2- 



so that 



a MME,2 



84c T s -M 

T\? 



(46) 



Relaxing the assumption of small QT/M, we can express the 
variance in a form similar to the CRLB ( fT4| > in Section |III-A| 

84<7 2 T S M 2 A „ iP" 



'MME, 2 



T 3 c 



1 + yiK- 



(47) 



where (1) K = QT/M, (2) higher j implies more accuracy 
and (2) jj are appropriate constants, e.g., first three terms can 
be derived as 71 = 24/35, 72 = 11/25 and 73 = 1581/6125. 
We skip the derivations for the sake of brevity. The variance 
can now be compared to the CRLB (given by ( fT4] i) and RLS 
(given by (|23|)). We can see that for small QT/M, expected 
variances of all of these techniques are proportional to the 
variance of the observation noise and inversely proportional 
to the cube of the number of samples. However, they have 
different multiplicative factors. The MLE has the best perfor- 
mance (factor of 48), followed by MME (factor of 84) and 
RLS (factor of 192). 

C. Notes on Implementation 

Some more remarks can be made about the performance 
of the MME estimator, given by ( p3] l, for a general value 
of m. For various moments m and small QT/M, we plot 
the multiplicative factors associated with a 2 T s M 2 /(T 3 c 2 ) in 
Figure |4] 

Choosing higher order moments decreases the variance 
initially before it starts increasing again. For example, for 
m = 3, the multiplicative factor is reduced to 64. The minima 
is attained for m ps 6 when the factor reduces to 53.4. Since 
MME is a sub-optimal method, its variance is always greater 
than CRLB which has a multiplicative factor of only 48. It may 
be noted that estimation using higher moments involves more 
computation and hence we have considered only m = 2, 3 in 
this paper for purposes of illustration. 

We remark on the computational efficiency of the esti- 
mator. Two components of our procedure require numerical 
computations - the integration of y t to obtain Y m x and the 
inverse of the function G m (r) to obtain r. Regarding the 
numerical integration, the trapezoidal method (e.g., [18|) is 
readily amenable for recursive implementation. Regarding the 
inverse, monotonicity and convexity of G m (r) (which can be 
established easily and skipped here for brevity) ensure that any 
iterative gradient based method can be used to compute the 
inverse. We have used Newton's Method in our simulations. 
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Fig. 4. Performance of MME based estimation in the transient region for 
different chosen moments. 



V. Illustrative Examples 
A. Experimental setup and online estimation 

Experiments were performed in the Building Energy and 
Environmental Systems Laboratory (BEESL) at Syracuse Uni- 
versity. We used a chamber of dimension 6.5 ft by 12 ft by 
10 ft high that was connected to an HVAC system with an 
independent roof top unit. The fresh air supply was set at a 
constant rate of Q = 28 CFM (cubic feet per minute). We 
use Q to know the ground truth but treat it as an unknown 
in our estimation procedure. To ensure well mixed conditions, 
a fan was used inside the chamber. CO2 concentration inside 
the room was measured (in parts per million volume) using 
Gray wolf sensor IQ-610, with a rated accuracy of ±3% read- 
ing and ±50 ppm. Constant inflow of CO2 gas was injected 
into the chamber using Alicat Scientific mass flow controller 
with a rated accuracy of ±0.8% reading and ±0.2% full scale, 
where full scale in this case was 20 SLPM (standard litres per 
minute). The source was a commercial gas cylinder containing 
100% C0 2 . The generation rate was maintained at 0.42 SLPM 
or 1.4 x 10~ 5 kg/sec, which corresponds to the rate at which 
a 170 pound person would produce CO2 while standing and 
performing a light task such as drawing |19|. This rate was 
also used as a benchmark in |8|. The sampling interval T s was 
chosen as 20 sec. The time constant of the growth curve can 
be estimated from these settings as M/Q w 83 samples or 28 
minutes. Injection of gas source into the chamber was started 
only 2 hours after closing the door, so as to ensure completely 
fresh air inside the chamber. The fresh air concentration Cq 
was taken as the mean of 30 samples before the source was 
injected and was measured to be 392 ppm. 

A sample online performance of the MME and RLS based 
estimators is shown in Figure [5] Since N = 1, the occupancy 
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estimator converges N — > 1 after sufficiently large number 
of samples. The top part of Figure [5] shows the actual noisy 
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observations. The steady state CO2 level can be computed as 
Cq+cN/Q w 922 ppm. As we can see, a reasonably accurate 
estimate can be obtained using fairly small number of samples, 
i.e., much before equilibrium is reached in the observed 
concentration. Also, the MME method provides estimates that 
appear less oscillatory and have lower variance compared to 
the RLS counterpart. We compare the statistical performances 
more rigorously in the next subsection by using Monte-Carlo 
simulations. 

B. Performance comparison using Monte-Carlo simulations 

In our experiments with the IQ-610 sensor, the steady 
state standard deviation was found to be approximately 8 
ppm. Hence for our simulations, we assume the measurement 
standard deviation = 10 ppm. To compare the performance 
of various schemes, we use the root mean square error of 

the occupancy estimate yVar(iV). We compare 4 schemes, 
namely MME with moments of order 2 and 3, RLS and 



MLE. The simulation setup is similar to that in Section V-A 



The results for a varying number of samples are shown in 
Figure [6] We observe that an MME scheme with the moment 
m = 3 performs better than the one with the moment 
tyi = 2. Both the MME schemes perform better than the 
RLS procedure. The relative performances of these schemes 
are partly predicted from the estimation variances derived in 
(fT4|> , (|23T> and (|47]i, and also the discussion in Section [TV-C| For 
quick reference, those metrics can be recalled as CRLB cx 
48a 2 /T 3 , a 2 MME3 cx 64a 2 /T 3 , a\ 1ME 2 cx 84cr 2 /T 3 and 
a\ LS cx 192a 2 /7* 

Next, we evaluate the accuracy of the transient-region 
expressions for variance derived in this paper. We use the 
third order expressions for the RLS (given by ( |23] l) and 
MME (given by ( |47] i) techniques. The constants ctj,f3j and 
7j for j < 3 were given earlier in Sections III-A III-B and 



|IV-B| respectively. We display both the theoretical performance 
predictions and Monte Carlo results in Figure [7] which are 
found to be extremely close. 
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Fig. 7. Theoretical approximations compared with Monte-Carlo performance, 
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C. Different metabolic rates 

An additional source of error in estimating occupancy from 
the CO2 generation rate is the fact that the volumetric rate of 
CO2 generation, c, depends on many factors such as Dubois 
body surface area Ad, the metabolic heat produced by the 
body Mjj and the respiration quotient RQ, which is the ratio 
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of CO2 exhaled to oxygen inhaled. Based on ASHRAE (e.g., 
0>0), the volumetric rate of CO2 generation per person can 
be written as 

_ 0.0028A D M H R Q 



0.23i? Q + 0.77 



(48) 



For an average size man Ap = 1.8 mt 2 . Metabolic rate 
Mjj depends strongly on the various types of activities and 
can range between 1 and 2 for the occupants of an office 
building (5). For an adult of average size engaged in light 
sedentary activity Rq is about 0.83. For our simulations, we 
have replicated an experimental classroom setup in [7], where 
the volume of the room is M = 6143 cu.ft., and the supply 
inflow -rate Q = 115 cfm. The sampling period was taken 
as T s = 2.5 min. We have considered population sizes of 5, 
10 and 20 and for all the cases we assume that each person 
generates CO2 with a metabolic rate uniformly distributed 
between 1 and 2. While estimating N, we assume that all 
the persons have a metabolic rate of M — 1.5. The results 
from 10 5 Monte Carlo trials are summarized in Figure [8] 
Because of the uncertainty in generation rates for each person, 
the estimates for the number of persons are inconsistent, i.e., 
do not converge to zero even for a large number of samples. 
The MME method is found to have better prediction capability 
for all cases. 
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continuous occupancy state as N', which is modeled as a 
random walk 



N[ +l =Nl + , 



with m l ~ d ' A/"(0,7 2 ), 



(49) 



where Ki are the disturbances that represent changing occu- 
pancy and 7 2 is the variance. The actual occupancy is modeled 
by the integer equivalent of N', i.e., 



Ni=[Nl\. 



(50) 



The carbon dioxide concentration is then calculated using 
the integral form of the mass balance equation ([3]) for each 
time step. The initial number of occupants is assumed to be 
Nq = 20 persons. Rest of the simulation setup is similar 
to that in Section [V-Cl i.e., M = 6143 cu.ft., Q = 115 
cfm, T s — 2.5 min. The estimation error is measured relative 
to the mean of the occupancy for the entire period and the 
results are displayed in Figure [9] for three different variances 
7 = 0.2,0.5 and 0.9. Since a random walk has a linearly 
increasing variance, the estimation error in Figure [9] increases 
with time. Also, the results show that the MME technique has 
better estimation capability than the RLS procedure. 
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D. Changing Occupancy 

One of our assumptions for both the RLS and MME 
estimation procedures has been that the occupancy is constant. 
In this subsection, we consider the case when the occupancy 
changes slightly during the duration of estimation. This is 
indicative of the situation when some attendees join or leave 
the meeting/ classroom. We consider a random walk model 
for occupancy, which was also used in JS). We denote the 



E. Timing and Sensor Reductions 

The improvement of estimation accuracy has direct implica- 
tions for other important system parameters, e.g., reduction in 
time to estimate and in number of sensors. From the transient 
region expressions of the variance, for small QT/M, we have 
for RLS ((23) and MME |47), 

a 2 RLS cx 192a 2 /T 3 , and a 2 MME 2 cx 84<7 2 /T 3 , (51) 

with identical constants of proportionality. In other words, for 
the same levels of acceptable variance, the RLS procedure 
would require -y/192/84 1.3 times the number of samples, 
i.e., 30% more than the MME procedure. To the ventilation 
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system, this will mean significant reduction in delay to re- 
spond. 

Also, it is well known that readings from multiple sensors 
can be fused to improve the overall accuracy ifTTI . In particular, 
N s independent sensors are known to reduce the variance of 
estimates by l/N s . If we consider a room instrumented with 
multiple and redundant sensors, then the RLS procedure would 
require 192/84 ss 2.28 times more number of sensors than the 
MLE technique. This may translate to significant savings in 
terms of instrumentation cost. 

VI. Conclusion 

In this paper, we have proposed a new approach for solving 
the problem of estimating the strength of a gaseous source in 
a room. We have assumed that the room is well mixed and a 
simplified mass balance equation is sufficient to describe the 
dynamics of the concentration over time. Also, it was assumed 
that the noise in the measurements is additive and Gaussian 
in nature. Since we are interested in making estimations as 
quickly as possible and much before equilibrium is reached, 
we have selected the initial section of the model dynamics 
as our operating region. We have performed a theoretical 
analysis of the estimation performances for our technique as 
well as another existing technique, and also obtained simplified 
expressions for this operating region. We have compared 
the performances of the two techniques using Monte-Carlo 
simulations and compared both with the lower bound. Our 
results clearly indicate the superiority of our approach. The 
proposed algorithm can potentially improve the performance 
of building control systems, since it can provide fast and 
accurate estimates of the strength of a gaseous source inside 
an indoor environment. 
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